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Abstract. The dynamics of two groups of molecular motors pulling in opposite 
directions on a rigid filament is studied theoretically. To this end we first consider 
the behavior of one set of motors pulling in a single direction against an external force 
using a new mean-field approach. Based on these results we analyze a similar setup 
with two sets of motors pulling in opposite directions in a tug-of-war in the presence 
of an external force. In both cases we find that the interplay of fluid friction and 
protein friction leads to a complex phase diagram where the force- velocity relations can 
exhibit regions of bistability and spontaneous symmetry breaking. Finally, motivated 
by recent work, we turn to the case of motility assay experiments where motors bound 
to a surface push on a bundle of filaments. We find that, depending on the absence 
or the presence of a bistability in the force-velocity curve at zero force, the bundle 
exhibits anomalous or biased diffusion on long-time and large-length scales. 
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1. Introduction 

Molecular motors are proteins which convert chemical energy into mechanical work. 
In many cases, relating to both in vivo and in vitro situations, they act together in 
large groups. Among the numerous examples are, myosin motors acting in muscles 
P^, kinesin motors pushing microtubules or myosin motors acting on actin filaments 
in motility assay experiments (see for example [2J E]), and the extraction of membrane 
nanotubes by kinesin motors [H EJ [6] . It is now well established that motors can exhibit 
a wide range of collective behaviors. Many times the collective behavior results in an 
oscillatory "like" motion [3, El El EE! EH] where the velocity changes abruptly between 
two distinct values [7J [12] . 

Frequently the setup is such that the motors all act together in a certain direction. 
This occurs, for example, in motility assays where the filaments have a well defined 
polarity. In some cases, however, the picture is different and two groups of motors pull in 
opposite directions. This is the case, for example, in muscle contraction, active gels with 
myosin minifilaments [T3l [TH [T5] . contractile ring that forms during cytokinesis, vesicles 
carried by both kinesin and dyenin along a microtubule bundles [16] and more. Recently, 
such a scenario has also been realized in motility assay experiments. In one set of 
experiments [17] a microtubule is acted upon by both NCD motors and homotetrameric 
kinesin-5 KLP61F in opposing directions. In another set, bundles of actin filament of 
opposing polarity are placed on a surface covered with myosin motors [18J. It is very 
common in such experiments to observe an oscillating like behavior where a velocity, 
say of a bundle of filaments, changes between two distinct values. 

So far theoretical studies of two classes of motors acting in opposite directions in 
a tug-of-war have focused on small groups of processive motors (which hardly detach 
from their track) acting on a fluid membrane [19]. In this paper we study theoretically a 
tug-of-war scenario for non-pro cessive motors acting in large groups on a rigid filament. 
We focus on motility assay experiments as discussed above, although many of our results 
can be easily extended to other scenarios. It is the central aim of this paper to analyze 
the different kinds of behavior that such a system can exhibit. To do this we generalize 
a discrete model, first introduced in [20], to incorporate two groups of motors pulling 
on a bundle in opposite directions. 

To this end we first revisit the usual scenario where a rigid filament is pulled by a 
specified number of motors of one type. We analyze the model through a new mean- 
field approach which allows a straightforward derivation of velocity-force relations, where 
the force is exerted by some external agent on the filament. Our mean-field approach 
gives rise to behaviors not observed in previous treatments of this system [201 EB 122] . 
We find four distinct types of force- velocity curves shown in figures [2(a) and [2(b). In 
particular the force-velocity relation exhibits distinct bistable behaviors which result 
from the viscosity of the fluid and the protein friction. The bistability manifests itself 
dynamically through an oscillatory like behavior where the velocity changes between two 
distinct values. Our analysis illustrates that there is a distinction between bistability 
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Figure 1: The setup we consider is of a bundle composed of many filaments set on 
a surface of motors. Note, that the analysis presented in the paper does not change 
if filaments of different polarity are mixed. In the inset a motor is shown with the 
transition rates between the attached and detached states. 



which arises due to a "fluid" viscosity studied in [21] from that of one which is caused 
by protein friction and studied in [23]. The model shows both. 

Building on these results we study the case of a bundle of filaments pulled in 
opposite directions by two groups of motors, with a specified number of each type, and 
acted upon by an external force (see figured] for a setup with no external force). Such 
a scenario could be realized using single molecule experiments. We find five different 
types of possible force- velocity relations with as much as four regions of bistability (see 
figures 11(a) and H](b)). Using existing data for myosin we discuss the force- velocity 
relation expected in a tug-of-war between two sets of myosin motors. We also present 
a systematic study of the dynamics of the system which result from the bistability. In 
the limit N — > oo the system can exhibit spontaneous symmetry breaking. 

We conclude by considering motility assay experiments. In contrast to the scenario 
discussed above the number of motors acting in each direction now fluctuates as a 
function of the location of the bundle on the substrate. Different behaviors are found if 
the system is in a bistable regime or not. When the system is not in a bistable regime, as 
the size of the bundle increases, its motion becomes irregular (see figure fT2l) . The bundle 
gets trapped for very long times at specific locations and in a large bundle limit (defined 
carefully in the text) the mean square displacement of the bundle grows as log A (t/r ). 
Here t denotes time and To sets a time scale. Moreover, at locations where the bundle 
remains trapped it displays an oscillatory like behavior despite the fact that there is 
no bistability. When the force-velocity relation is bistable the motion is also irregular. 
However, now in the large bundle limit the mean square displacement of the bundle is 
linear in time (diffusive) and the bundle is expected to have a mean, non-zero, velocity. 
Our results provide an explanation for the experiments of [18] without resorting to their 
assumption of cooperative unbinding of the motors. 
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The structure of the paper is as follows. In section [2] we define the model. In 
section [3] the model is analyzed using a mean-field approximation. In section H] the 
dynamics of the motion are examined. Section [5] studies the case where the number of 
motors varies randomly in space. Finally we summarize in section [6j 

2. The model 

The scenario of interest is of two sets of motors pulling in opposite directions on a 
common backbone - a tug-of-war. One set pulls in the "plus" direction with a force F + 
while the other pulls with a force F_ in the "minus" direction (see figured]). The motion 
at low Reynolds number is governed by the force balance equation where the force on 
the filament is countered by the viscous drag of the filament, namely 

rjv = F + + F_ + F ext . (1) 

On the left hand side of the equation rj is the viscosity and v is the velocity of the 
filament. On the right hand side are the forces acting on the filament and we include 
a possible contribution from an external force, F ext . We assume a completely rigid 
filament. To find F± a microscopic model of the motors is needed. 

The model we use is based on the one introduced in [20]. Illustrated in figure [U it 
consists of two motor states denoted by a and d. In the a state the motor is attached 
to the filament, while in the d state it is detached. The transition rate between states 
d and a (a and d) is denoted by k a (kd). As shown in figure (TJ we assume that when 
the motor binds to the filament it enters a tense state, where some "spring-like" degree 
of freedom, x, is extended from an equilibrium position x = 0. The motor then exerts 
a force as the spring relaxes and finally completes the cycle by detaching from the 
filament. In general we expect kd to increase with the tension on the motor. Following 
|20j . we take kd = uj^exp (a\x\) which is consistent with the usual Kramers form of 
rates with a = Kl/ksT. Here if is a spring constant, I is a microscopic length, T is 
the temperature and ks is the Boltzmann constant. In addition to being physically 
motivated, as stressed in [20] , simpler forms, for example monotonically decreasing with 
x, do not yield an oscillating like behavior which is the focus here. 

To analyze the model we employ a mean-field approximation. This is done by 
writing self-consistent equations for the force generated by the motors and the transition 
rates. To this end we relate the displacement x (of the "spring-like" degree of freedom) 
to the velocity of the filament, v. We expect the approach to hold for a large number 
of motors, when the fluctuations in the velocity are negligible. The self-consistent 
equations are obtained as follows. We denote the velocity of the filament by v, the 
time since the motor attached to the filament by t and the initial extension of the spring 
after attaching to the filament by xq. The rate kd(x) can be expressed using xq, v and 
t, namely k d (x 0l v,t) = kd(x) = kd(x — vt). The probability density of detaching at 
time t is then p (t, v, x ) = exp ( — f* kd (x , v, t') dt' ) k d (x , v, t) so that a self-consistent 
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(2) 



Here q(xo) is the probability density of the springs attaching to the filament with an 
extension xq. Note that within our approximation we replace the several time scales 
present in p(t, v, xq) by a single one. The self-consistent rates can be easily obtained for 
different choices of q(xo) numerically. 

We find that for generic choices 1/kdiv) is peaked around some value v — v and 
decays monotonically to zero on both sides of the peak. Fast velocities (positive or 
negative) imply a fast growth of kd(x). This in turn makes kd(v) larger. The peak 
around the finite value of v results from the finite positive average of q(xo) which is 
assumed. Moreover, we find that the exact form of the function is unimportant for the 
qualitative results expected from such a simple model [24| (this will become evident 
later). To this end, we use a simplified analytical form, kd (v) — k ((?; — v ) 2 /W + l), 
which captures all the important features of the exact form derived by the procedure 
above. Here k is the binding rate at v — Vo and W is a scale parameter with units 
of velocity squared. A similar procedure is used to define k a . Since in the detached 
state the spring is not stretched, it is easy to see that the above procedure leads to an 
attachment rate which is independent of v (see also [201 122] ) . 

Next, we need to specify the force exerted by the motor. This force, of course, 
varies with time. Within our approach we replace the time dependent force by its 
average, (f(v)), over the attachment time to the filament l/kd(v). While this can be 
done formally it is easy to see that to linear order in v, (f(v)) = G — ^v. Positive 
velocities tend to initially release the tension in the spring decreasing the force. The 
term jv is the leading order behavior of the protein friction [22] which arises due to 
the energy dissipated when a tense motor unbinds. To see this, consider motors that 
are unable to actively generate force, so that G = 0. Clearly, — is the force resulting 
from the elastic element of the bound motors being stretched by an external force. In 
general 7 can have a non trivial, details sensitive, dependence on the velocity. However 
for large velocities motors are quickly detached from the filament and the effect of the 
friction becomes less important. This will be shown in the treatment of the collective 
behavior of the motors presented below. 

Before turning to the two sets of motors problem it will be useful to first analyze 
the single set problem, F_ = 0, as our solution to the two set problem relies on it. 

3. Mean-field analysis 

3.1. One set of motors. 



In this section we will examine the case where a single set of motors operates against 
an external force, namely F_ = 0. In principle, since the model constitutes a one-step 
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process over the number of motors in the a state, N a , a formal solution may be written 
for the steady state. However, it is more instructive to examine a mean-field solution 
that is easily generalized to two sets of motors. Based on the rates defined above it is 
straightforward to write mean- field equations for Pj = Ni/N, the fraction of motors in 
each state 



Here iVj is the number of motors in state i = d,a and iV is the total number of motors. 
Note that a single motor is coupled to the rest of the motors through the rate kd(v) 
which has a non-trivial dependence on the velocity, v. The mean-field approximation, 
on top of the approximations described above, neglects correlations between kj(v) and 
P a . Furthermore, since the motors are identical the same equation holds for all motors. 

The stationary solution of these equations is easily solved and along with the 
normalization condition Pd + P a = 1 yields 



Note that the expression depends on v, the velocity of the filament. To obtain v as 
a function of, say f ex t = F ex t/N, one then uses the solution self-consistently in the 
force balance equation, setting F + = N a (f(v)). To a first approximation we take fj 
proportional to the length of the filament. Assuming the motors are evenly spaced out 
the overall viscosity can be expressed using the total number of motors N and rj, the 
viscosity per unit distance between motors, such that fj = Nrj. This results in: 



with P a {v) given in (j3J). The right hand side of the equation is a sum of three terms. The 
first, due to the viscosity of the filament, is monotonic in v. The second results from the 
protein friction and the third from the force exerted by the motors. The last two terms 
multiply P a {v) which due to the functional form of kd{v) is a non-monotonic function of 
v. As we show this implies that in some ranges of the parameters there is a region where 
for every value of f ext there are three solutions for f ext (v). Following [12] we take this 
as evidence for bistability. In fact, we demonstrate below that the competition between 
the three terms, with two non-monotonic, can lead to a rather rich behavior with four 
distinct phases. For now we focus on the steady state solutions. Later on we analyze 
the dynamics which arise in the bistable regimes. 

3.2. Phase diagram 

In this section we classify the different possible force-velocity curves. Asymptotically, 
for large f ext most of the motors are detached and the velocity is roughly v ~ fext/V- 
For smaller forces, depending on the parameter values, the contribution from the non- 
monotonic terms can be important. We find four distinct regimes which qualitatively 
depend on 7, G and r\ as follows: 



d t P a = P d k a - k d (v)P a . 



(3) 




(4) 



fext (v) =T)V+ (jV - G) P a (v), 



(5) 
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(i) Large r], small 7 and small G - No bistable regime: Here the first term in §5$) 
dominates due to the large viscosity so that the velocity changes monotonically with 
the external force. The stall force, defined as v (f ex t) = occurs for f ext < 0, or in other 
words the velocity is positive when there is no external force acting on the motors (see 
figure [2]( a)). 

(ii) Small r\, small 7 and large G - Single bistable regime at f ext < 0: Here the third 
term, due to the force exerted by the motors, is large enough so that its non-monotonic 
behavior becomes important. This is seen in figure EJ^b) where the force-velocity curve 
shows a region with three possible values of v for a given value of f ext . In particular 
there is a region where d v f ext (v) < for one solution and d v f ext (v) > for the two other 
solutions. Since d v f ext (v) < implies a negative mobility this suggests that the solution 
is unstable. The difference between the two stable solutions is manifested through the 
number of motors in state a. The solution with the larger velocity has more motors 
in the attached state enabling it to counter the external force. The solution with the 
smaller velocity has most of the motors detached, resulting in a negative velocity. This 
can be seen by setting 7 = in (jSJ) and is illustrated in figure H](d) where we plot 

Pa (fext) • 

(Hi) Small G, small r\ and large 7 - Single bistable regime at f ext > 0: Here only the 
second non-monotonic contribution due to protein friction induces a bistable regime. 
This occurs as long as rj is not too small that the non-monotonic behavior due to the 
force exerted by the motors becomes unimportant. Now the force- velocity curve has 
a single bistable region located at f ext > while the velocity varies continuously for 
fext < (see figure [2](a)). Unlike (ii) both stable solutions have a positive velocity. 
(iv) Small 7], large 7 and large G - Two bistable regimes: Here both non-monotonic 
terms, due to the protein friction and the force exerted by the motors become important. 
Interestingly, this can lead to two distinct regions of bistability. In contrast to the 
bistable region due to the force exerted by the motors (ii), in the additional bistable 
region both of the stable solutions have a positive velocity. This is illustrated in 
figure [21(b) where it is seen that increasing 7 adds a second bistable regime located 
at positive forces. 

For the specific choice we make here for kd(v) the results can be understood as 
follows. Due to the symmetric choice of the rate kd(v), the fraction of attached motors 
has the property P a (v + v ) = P a (—v + v ). Therefore, the force of the motors can be 
expressed in terms of symmetric and anti-symmetric functions relative to vq, namely 

(G - iv)P a (v) = P a (v)(G - v ol ) - 1 P a (v)(v - vO). (6) 

The first term on the right hand side is symmetric with a single extrema and therefore 
contributes a single non-monotonic region. The second term on the right hand side is 
anti-symmetric and has two extrema located symmetrically around Vq contributing two 
bistable regions. Note that when G = Voj there is only an antisymmetric contribution. 
We stress that while this argument is specific for our choice of ka(v ) the general structure 
is unchanged for other choices of non- symmetric ks,(v). 
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Figure 2: Results for a single set of motors, (a) force- velocity curve of types (i) and 
(iii). Here 77 = O.lpNsec/ fim. (b) Force-velocity curves of types (ii) and (iv). Here 
77 = O.OlpNsec/ jjtm. (c) Bistable ranges of f ext as a function of 77 (shown in blue and 
red), (d) Fraction of attached motors with parameters as in (b). In (a) and (b) the 
dotted line signifies the regions of the solutions where d v f ex t(v) < 0. In all figures 
k a /K — 1/5, G = 5pN, W = 10(fim/sec) 2 , v = 2fim/sec. 7 is given in the graph in 
units of pNsec/ fim 



The different parameters can in principle be controlled, to some extent, 
experimentally to observe possible transitions between the different regimes. For 
example, 77 can be controlled by the density of motors along the filament and 7 can be 
controlled to some extent by the ATP concentration. It is well known that increasing 
ATP concentration increases the unbinding rate [25] . Therefore large values of 7 would 
correspond to small unbinding rates of the motors and hence small ATP concentration. 
G might be tuned by changing the neck region of the motor. In figure [2](c) we illustrate 
two transitions controlled by the viscosity 77 — one from two bistable regimes to a non- 
bistable regime and one from a single bistable regime to a non-bistable regime. The 
filled areas represent ranges of f ex t for which there is bistable behavior for the first 
transition (blue) and the second transition (red). For every value of 77, the region is 
calculated from the two adjacent extrema of the mean-field f ex t (v). For large enough 
viscosity 77, the velocity v changes continuously with the external force, f ext . As 77 is 
decreased beyond a threshold value a bistable region emerges and the range of force 
values it encompasses grows. 

It is interesting that the simple model presented above accounts well for the 
measured force- velocity curves for myosin II. In figure [3] we show a fit of the mean-field 
solution to the data of reference [11] along with numerical simulations of the model. 
In contrast to the mean-field solution, the experimental and numerical data show that 
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Figure 3: A fit of the theory to experimental data. Solid dots (•) are data extracted 
from reference [11]. The solid line is the results of the mean-field theory and empty 
circles (o) represent the average velocity obtained from the simulations of the model 
using the fitted parameters of the mean-field results. The dashed line represents the 
predictions for tug-of-war between two equally sized groups of myosin-motors. 

the positive velocity branch is stable up to a certain force. In other words, up to this 
force the motion of the filament with negative velocity is hardly observed. We shall 
return to this point. It results from of the dynamical aspects which are not captured 
by the mean-field solution that will be discussed in a separate section and appears in 
our numerical simulations of the model (shown in the figure). To fit the data we used 
n/k a = 1/10, G = 4.5pN, vo = 1.5/im/sec, 7 = 0.28pNsec/ fim, W = 16(/im / sec) 2 
and r) = O.OOGpNsec/ fim, which are in the range of values measured independently in 
references [25] and [26] (see also references within). The experiments were performed 
using an external voltage and we rely on their conversion between voltage and force 
per motor. While in [11] the points with negative velocity are quoted to be stable our 
model predicts them to result from the motion of the filament in two opposite directions. 
Clearly more data is needed in this region to improve the fit. Moreover the fit can be 
improved by using more elaborate forms for kd(v). 

The above analysis interpolated between two regimes which were previously studied. 
In one class of models a joint rod and protein friction term was accounted for in the 
form of the first term on the right-hand side of (13) [2TJ. In another only protein friction 
in the form of the second term on the right-hand side of flS]) was considered [20]. We 
show that both lead to distinct bistable regimes which can lead to a richer behavior than 
previously discussed. Our focus in the paper, however, is a tug of war scenario. As we 
illustrate, using the approach developed above, the analysis becomes straightforward. 
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3.3. Two sets of motor engaged in a tug of war 

We now turn to discuss the possible force-velocity curves for two sets of motors pulling 
one against another in a tug-of-war. Again we focus on the steady-state mean-field 
predictions and later discuss dynamical aspects. The setup we consider is as follows. 
Two (or several) connected actin filaments of opposite polarity form a bundle which 
is set on a substrate of motors. In this section we assume that the number of motors 
which operate on the bundle is constant. Thus, the filament (or filaments) with positive 
polarity is acted upon by N + motors pulling in the positive direction while the filament 
(or filaments) with negative polarity is acted upon by N~ motors pulling in the negative 
direction (see figured]). 

Motivated by the experiments of reference [IS], and for simplicity, we assume that 
both sets of motors, pulling in opposite directions, are identical. The definition of the 
rates of the model then follows as above but now with a modified force balance equation 
which accounts for both sets of motors. Defining iV+ (N~) to be the number of attached 
motors pulling in the positive (negative) direction the force balance equation now reads: 



Here N + + N~ is the total number of motors, iV+ (G — jv) is the force produced by the 
motors pulling in the positive direction and — N~ {G + jv) is the force produced by the 
motors pulling in the negative direction. Finally F ext is an external force (the overline 
is used to distinguish this from the case of a single set of motors). As in the previous 
section the total viscosity is taken to be proportional to the total number of motors. 

Again we treat these equations within a mean-field approach replacing N^ 1 by 
N^P^ 1 . Here P^[ (P^) is the fraction of motors in the attached state pulling in the 
positive (negative) direction. Clearly, these satisfy (J4]) so that (v) and P~ (v) are 
related through P+ (v) = P~ (—v). Using this with the expression for f ex t(v) defined in 
(JSJ) the force-velocity relation takes the form: 



Note that the relation only relies on the anti-symmetric nature of the force exerted by 
the motors and additivity of the viscous terms. In the following we consider the case 
N + = N~ = N, for which F ext (v = 0) = 0. For this case the force-velocity curves, in 
terms of the normalized external force f ex t(v) = F ext (v)/N, are independent N. The 
force- velocity curves can be easily obtained using (JSJ by creating antisymmetric combi- 
nations of the curves obtained for one set of motors. 

3.4- Phase diagram 

By following the above procedure we find four distinct generic structures for the force- 
velocity curve. Similar to the one motor case these arise from an interplay of the various 
viscous terms and the force exerted by the motors. In contrast to the single motor case 
the generic curves cannot be easily classified according to the values of the viscosities 



V v (N + + N~) = N+ (G - jv) - N~ (G + jv) + F ext . 



(7) 



F ext (v) = N + f ext (v) - N~f ext (-v) . 



(8) 
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and the force exerted by the motors. 

(I) Monotonic force-velocity curve: Here the force- velocity curve does not have any 
bistable regions (see figure H](a)). Such behavior can arise under several situations. The 
simplest one involves an antisymmetric combination of a monotonic force- velocity curves 
for a single motor (regime (i) of the previous section). However under certain conditions 
it might also occur by antisymmetric combinations of any of the other regimes. This 
will occur when the non monotonic regime of one curve is weak enough so that when a 
monotonic contribution is added to it, it becomes monotonic. In all these cases, in the 
absence of an external force the average velocity of the rod is zero. 

(II) Single oscillating regime, centered around f ex t = 0: This can occur only when we 
combine single motor curves from regime (ii) or regime (iv) of the previous section. 
The latter combination requires that the non monotonic behavior due to the protein 
friction is canceled by the antisymmetric combination. In this case at zero external 
force the system exhibits two stable solutions with velocities equal in magnitude but with 
opposite sign (see figure IH(a)). For the solution with positive velocity > N~, while 
for solutions with negative velocity iV+ < N~ (see figure H](d)). Note that similar force- 
velocity curves have also been shown to exist theoretically for completely symmetric 
motors [12] . 

(III) Two non-monotonic regimes This occurs when single motor force-velocity curves 
from regime (ii), regime (Hi) or regime (iv) of the previous section are combined. Now 
the force- velocity curve has two symmetric bistable regions at positive and negative 
forces (see figure H^b)). In the absence of any external force there is a single solution 
with zero velocity. In this case the fraction of attached motors is equal and relatively 
large. 

(IV) Three bistable regimes: This can occur only by combining force-velocity curves 
from regime (iv) of the previous section. Now there is one bistable regime around zero 
external force and two antisymmetric bistable regimes at positive and negative external 
forces. 

(V) Four bistable regimes: As evident from the above constructions it is also possible 
to find regimes where there are four bistable regimes by suitable combinations of curves 
for one set of motors. Two are at small positive and negative values of the forces and two 
are at large ones. In fact, in some cases, the two small force regimes may be located one 
on top of the other. This leads to a regime around zero force with four stable velocities. 
By explicitly plotting the curves we find that both of these require a careful fine tuning 
of parameters. We therefore do not expect them to appear under generic conditions. 

The transitions into the bistable regimes as a function of viscosity are illustrated in 
figure H^c) for various parameters. Again we see that increasing the viscosity smooths 
out the force- velocity curve until it becomes completely monotonic. In addition P^(fext) 
is seen to have a complex structure, as shown in figure H^d), with the same number of 
bistable regions as the force-velocity curve. 

In the case of two equally sized groups of motors, a simple criterion for bistability 
at fext = can be found from the results for a force-velocity curve of a single set of 
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Figure 4: Results for two equally sized sets of motors pulling in a tug of war. (a) 
Force- velocity curves of types (I) and (II). Here rj = 0.01. (b) force-velocity curve of 
types (III) and (IV). Here rj = 0.04. (c) The range of values of the external force where 
there is a bistability as a function of rj (marked by green, blue and red areas), (d) The 
fraction of attached motors pulling in the "positive" direction as a function of external 
force, P + (fext)- In all figures k a / n = 1/5, G = 5pN, W = 10(fim/sec) 2 . 7 is given in 
the graph in units of pNsec/ fim. 77 is specified in the graph in units of pNsec/ fim. For 
the green curves (I), (III) v = l/im/sec, for the red curves (II) v = 0.6[im/sec and for 
the blue (IV) curves vq = 4fim/sec. 



motors. This is done by requiring d v f ext \ v =o < 0, which can be rewritten using (JSj) : 

d v fext\v=o = 2d v f ext (v)\ v=0 < 0. (9) 

Namely, the slope of the single motor set force-velocity curve at stall force has to be 
negative. On a single set of motors this could be measured from the slope of the force in 
constant velocity experiments [23J around zero velocity. Constant velocity experiments 
are required as the zero velocity solution might be unstable. 

We now turn to discuss the case where N + ^ N~ . The force- velocity curves are once 
again found using the antisymmetric combinations of the single force- velocity curves (see 
equation (jSJ)). When this is done the force- velocity curve is no longer antisymmetric. 
Changing the ratio N + /N~ leads to a continuous change of the force- velocity curve as 
shown in figure [5j As expected, as the ratio of the number of motors increases the curves 
change continuously from the behavior of two sets of motors to that of a single set of 
motors. In general this leads to a loss of possible bistable regions of the force-velocity 
curve. 

Finally, we note that using the parameters from the fit of figure [3] it is possible to 
predict the force-velocity curve for two groups of myosin motors. We expect a force- 
velocity curve of type (//), namely at zero force there are two possible velocities. Note, 
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fext \PN] 

Figure 5: Results for different numbers of motors of each type. Shown is the evolution 
of the force- velocity curve as the ratio N~ /N + (given in the legend) changes. f ext is 
normalized by iV + . Here n/k a = 5, G = 5pN, W = 10 (/xm / sec) 2 , v = 3fim/sec, 
7 = l/im/sec and rj = 0.02/im/sec 



however, that care must be taken with this conclusion since we applied a specific (simple) 
model to fit the experimental data. For a clear conclusion the procedure explained above 
should be performed on data using constant velocity experiments (for example, using 
single molecule techniques). 

4. Dynamics 

We now discuss the dynamics of two sets of motors pulling oppositely in a tug-of-war. 
As above we assume a setup where N + motors pull against iV - motors. We focus 
on the generic case (as we argue below) that N + m N~ and f ext is small. The more 
conventional setup of disordered motility assays will be discussed in the next section. 

The mean-field treatment of above does not account for the dynamics, and only 
provides possible steady state solutions. To account for the dynamics we expand the 
master equation of the process in powers of l/N^ into a Fokker Planck equation. These 
results, and all other numerical results presented in the paper, are verified using Monte- 



Carlo simulations, described in Appendix A 



As we show below two generic behaviors which depend on the presence of a 
bistability in the mean-field solution are found. When there is no bistability (regimes 
(I) and (HI)), as intuitively clear, the motion on long time-scales is a biased diffusion. 
The bias vanishes for N + = N~ at f ext = 0, while for N + ^ N~ the bias vanishes for a 
non zero value of f ex t- Note however, that on short time scales the motion may display 
an oscillatory like behavior as seen in figure ITTT b). The velocity distribution function 
can not be fitted using a single Gaussian function (see in figure HIT a)). 

When a bistability is present (regimes (II) and (IV)) the behavior is more 
interesting. We find stochastic transitions between the two mean-field solutions. The 
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motion on each of the two mean-field solutions is a biased diffusion with the bias dictated 
by the corresponding velocity. Each is characterized by an average dwell time, r + and t_ 
for positive and negative velocities respectively. These are defined by the average time 
in which one solution changes to the other. Of course on long time scales the behavior is 
still a biased diffusion. The analysis below shows that the average dwell times increase 
exponentially with the number of motors. The prefactor in the exponential depends on 
the exact value of the external force (This is very similar to the exponential time scales 
found for the symmetric motors in [12]). For N + = N~ = N the ratio of the dwell times, 
7+/t_, is exponentially large (small) in the number of motors for f ext > (f ext < 0), 
namely t+/t_ ~ exp(6N f ext ) where 9 is a positive constant. This implies that for a 
large number of motors the transition between the two branches becomes very sharp. 
Specifically, in the N — > oo limit at f ext = the system spontaneously breaks between 
the two directions of motion. As discussed in [12] a similar picture also holds for a 
single set of motors around the force regime exhibiting bistability. Note, that when no 
bistable region is present the corresponding dwell times of the oscillations like behavior 
(see figure [TTTb)) has a weak dependence on the number of motors. 

To derive these results we start by analyzing the case of a single set of motors 
which is then easily generalized to two set of motors. Our model for one set of motors 
is a one step process with N distinct states defined by the number of attached motors, 
N a . The transition rate from a state with N a motors to one with N a + 1, is given by 
9Na = k a (N — N a ). The transition from N a to N a — 1 is given by r Na = N a k d {v). The 
master equation then reads 

d t p(N a ) = -(r Na + g Na )p(N a ) + r Na+lP (N a+ i) + 9 N a -A N a-i) , (10) 

where p{N a ) is the probability of having N a motors attached. Next, we define q = N a /N 
so that for every N, q takes values in the range to 1. For brevity we use the operator 
E, defined through its operation on a function f(q), namely E/(g) = f(q + l/N) and 
•/(<?) — /(? ~ V-W)- It i s eas y t° see that the master equation using these notations 
is given by 

d t p(q) = (E - l)p(q)r(q) + (E' 1 - l)p(q)g(q), (11) 

with r(q) = r(Nq) and g(q) = g(Nq). Next the operator E is expanded in powers 
of l/N [2^ keeping terms up to second order so that: E = 1 + j^d q + jw^^iq an< ^ 
E _1 = 1 — jfd q + d qq . This gives the Fokker-Planck equation: 

dtp(q) = -d q F q p(q) + ±^d qq (D g p(q)). (12) 

The right hand side of the equation has two terms: a drift term defined through, 

F q = ^(g(q) - r(q)) = k a (l - q) - qk d (v), (13) 

which for a constant f ext is independent of N. Note that requiring F q = yields the 
steady-state mean-field equations. Therefore, the previous mean-field solutions serve as 
extrema of the effective potential of the Fokker-Planck equation. It is easy to verify that 
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solutions which satisfy d v f ext < correspond to maxima while those with d v f ext > 
correspond to minima. The second term is a diffusive term, with 

D q = jj(gNq + r Nq ) = k a (l -q) + qk d (v), (14) 

independent of N for a constant f ext . In both the expressions of F q and D q the velocity 
is given by the force balance equation: ryv = q{G — jv) + f ext . Note the overall 1/N 
term multiplying the diffusive term. A standard Kramers analysis implies that the time 
needed to cross the barrier from one minima to the other is given by r ~ exp(aiV) 
where a is a positive constant which depends on the functional form of D q and F q . The 
analysis shows, similar to [12], that the time increases exponentially with N. 

Next we extend the above analysis to two sets of motors. We focus on the case 
fext = 0. Repeating the same procedure used to obtain ( Tl2l) but now with the two 
coordinates q = iV+ /N + and s = N~ /N~, one obtains: 

d t p = - d q (F q p) - d q (F s p) 

+ ^+d„ (D q p) + ^-d ss (D s p) . (15) 

F q and D q are defined as above and F s and D s are given by: 

F s = k a (l - q) - qk d (-v), 

D s = k a (l-q) + qk d (-v). (16) 

On the right hand side of (|15p there are now two drift terms and two diffusive terms, 
one for each direction. This equation is augmented by t](N + + N~) = qN + (G — ^v) — 
sN~(G + 71'). To analyze the equations we consider first the case N + = N~ . In this 
case it easy to see that, due to the symmetry of the force balance equation, F q , F s , 
D q and D s become independent of A r± . As in the case of a single a set of motors, the 
steady-state mean-field solution corresponds to F q = and F s = 0. 

The effective force-field may have either one or two stable stationary points 
that correspond to the number of the stable mean-field solutions (we ignore the rare 
possibility of four bistable regimes). We begin by discussing the latter and argue that it 
leads to "oscillations" between the two solutions. In figure IHa] the effective "force-field" 
caused by F q and F s for general values of s and q is shown. One can clearly see two 
stable minima and a single stationary unstable point (saddle point). The trajectories 
dtq = F q (q, s) and dts = F s (q, s) are exact in the limit N ± — > oo. It is straightforward 
to argue that r + = t_ = exp(/5A r+ ) = exp(/9iV~). Here (3 is a positive constant and r + 
and r_ are the typical dwells time in the minima corresponding to the motors pulling 
to the right and left respectively. 

Generically iV + and iV~ are not equal. In particular one expects for a typical 
bundle, composed of a random assortment of filaments with opposite polarity (see 
figured]), N + — N~ ~ \/N + + N~ . In this limit the mean- field solutions are modified 
from the case N + = N~ by terms of the order 1 / ^/N + + N~ , vanishing in the large A r± 
limit. In contrast, we find that the dwell times are very sensitive to the difference in 
the number of motors of each type. In particular we find r + /r_ ~ exp(/3A r+ — /3N~). 
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Figure 6: The effective force-field that governs the dynamics when there is a bistable 
region, (a) The effective force-field in the s , q plane plotted when N + = N~ . The 
three stationary points are marked by circles. Two of them are stable points while the 
third is a saddle point, (b) The additional contribution of the force-field that result for 
N + > N~ . Here we set: n/k a = 5, W = 10 (fim/sec) 2 , vo = 2 fim/sec, G = 5 pN, 
7] = 0.02 pNsec/ fim and 7 = 0.4 pNsec/ fim. 



Namely, when oscillations are present a small relative difference in the number of motors 
can lead to a strong asymmetry in the dwell times. 

To see this dependence we consider the case A r± = iV(l ± e) with e of the order 
1 / \J~N . The force-balance equation now gives v an explicit dependence on the values of 
iV + and AL. This in turn leads to corrections of the order of e to both F q and F s , which 
are easily found to first order: 

5F q = - e2 g— - v ) 5v, 

SF S = - 6 2s^-(v + v )5v, (17) 

Gq + Gs (q- s) 2 G 7 

dv = 7) . 

2rj + 7g + 7s (2rj + 7 g + 7s) 

The contributions of these terms, as seen in figure l6bl increase the transition 
probability from one minimum to the second and decrease the transition in the opposite 
direction (depending on the sign of e). The analysis follows as above and we find 
t + oc exp(aiV + neN) and r_ oc exp(aiV — x e N) with a, k and x positive constants. 
This implies that 

r + /r~ ~exp((/€ + X )eiV), (18) 

so that the ratio of the dwell times increases exponentially with the difference between 
the number of motors of each type. Note that the analysis relies on a relative small 
difference in the number of motors and in general holds as long as e is small. When 
this is not the case the results might change. In particular, the mean-field solutions, 
as discussed in the previous section, might not show a bistable behavior. To verify the 
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Figure 7: The ratio of the dwell times, r + /r_, in a bistable region when N + ^ N~ and 
there is no external force as obtained from numerical simulations. As can be seen to a 
good approximation r + /r_ oc exp(aN + — aN~), where a is a positive constant. Here 
n/k a = 1, W = 10 (fim/sec) 2 , v = 3 fim/sec, G = 5 pN, r\ = 0.0125 pNsec/ \im and 
7 = 0.125 pNsec/ \im. 

predicted behavior we have carried out numerical simulations. In figure [7] we plot the 
ratio t + /t_ and show that indeed it behaves as discussed above. 

A similar effect occurs when iV + = N~ = N and an external force is added. To 
first order in f ext the effective forces F q (s, q) and F s (s, q) have the additional terms: 



It is straightforward to argue that t + /t_ oc exp(6Nf ext ), where 9 is positive constant. 
This is seen in figures E] and O Note that this implies that for large A r± the velocity, 
v(f ext ), switches sharply between the positive and negative mean-field branches (see 
figure [3J. In the limit N — > oo at f ext = the system exhibits spontaneous breaking of 
the symmetry between the two directions of motion. 

We now turn to discuss the case where there is no bistable region and N + = N~ . 
Since there is only one steady-state mean-field solution, the effective force-field has only 
one stationary point. Clearly, the motion on long time scales is therefore diffusive. 
However, on short time scales the motion may display a weak oscillatory like behavior. 
Figure [10] shows an effective force that leads to this behavior. At the vicinity of the 
stationary point the force on the constant q + s line is smaller than the force on the 
q = s direction. This induces anticorrelations between q and s which increases the 
probability to find a non zero velocity. Typical traces of the displacement of the bundle, 
x(t), appear almost bidirectional and the velocity distribution deviates from a Gaussian 
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Figure 8: The ratio of dwell times, r + /r_, in a bistable region as a function of N + + N~ 
for N = N + = N~ and a constant f ext as obtained from numerical simulations. 
As shown to a good approximation, r + /r_ oc exp(0f ext N). Here k a / k — 1, W — 
10 (fim/sec) 2 , v = 2.5 fim/sec, G = 5 pN, r\ = 0.01 pNsec/ fim, 7 = 0.5 pNsec/ \xm 
and f^t = 0.12 pJV. 




0.05 0.1 0.15 0.2 



fext \pN] 

Figure 9: The ratio of dwell times in a bistable region for a constant iV = iV + = iV - 
as a function of f ext as obtained from numerical simulations. As can be seen, to a good 
approximation, r + /r_ oc exp(0f ext N). All the parameters are the same as in figure M 
except for vq = 3fim/sec. 



distribution as seen in figure HIT a). Note that for other choices of parameters, q and s 
can be positively correlated, yielding a Gaussian distribution. Finally, as clearly evident 
in the case N + 7^ the average velocity is no longer zero. 

In summary, irrespective of the presence of the bistability in the mean-field solution 
for a fixed and A^ - , the motion on long-times and large-length scales is a biased 
diffusion. When a bistability is present the motion will exhibit oscillations in the sense 
described above and the crossover to the final biased diffusive behavior is expected to 
occur on time scales which are exponentially large in the number of motors. 
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Figure 10: The effective force-field when there is a single steady-state mean-field solution 
plotted in the q, s plane. Note that there is only one stationary stable point. The shape 
of the effective force-field induces anticorrelation between q and s. The parameters are 
the same as in figure l6al except for 77 = 0.01 pNsec/ \im and 7 = 2.5 pNsec/ \im. 
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Figure 11: Dynamics in the absence of a bistability that result in an oscillatory like 
motion, (a) The velocity probability distribution obtained from numerical simulations, 
(b) A representative trajectory x(t) showing an oscillatory like motion. The parameters 
used are the same as in figure [103 



5. Disorder 

So far the number of motors pulling in each direction A r± was taken to be a constant 
independent of the displacement of the filament. In typical motility assays the filament 
is moved around by motors which are bound to a substrate. Since the motor distribution 
is expected to be random this implies that in general A r± depends on the location of the 
center of mass of the filament, its orientation and possibly time. In a typical experiment 
one would expect A r± to be a random variable. This results, for example, from, as 
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stated above, inhomogeneities of the motor density, randomly orientated motors and 
the motion of the motor head. In our analysis we ignore orientational changes of the 
bundle and assume that the motion is along a single axis. In motility assay experiments 
one dimensional motion is expected to occur on short time scales or by restricting the 
motion by methods used in [TT| . In fact some experiments show a motion very close to 
one dimensional [T8] . 
Two extreme cases may be considered: 

(a) N^ 1 is only time dependent. This could result, for example, from extremely flexible 
motors and a bundle composed of a fine mesh of filaments. This allows each motor 
to bind to a positive and negative filament with equal probability. The motion in this 
case on long-time scales and large-length scales will clearly be biased diffusion with the 
average velocity determined by the average of N + — N~ . A non-zero average N + — N~ 
could result, for example, from a difference in the amount of filaments of each polarity 
which compose a bundle. To conclude, in this case the additional randomness does not 
change the qualitative motion however it could change the diffusion constant and the 
bias. 

(b) N ± is only x dependent. Here x denotes the location of the center of mass of the 
bundle. This can occur, for example, if each motor can bind only to one of the directions 
- either the plus oriented filaments or the negative oriented filaments and will naturally 
be the case when the bundle of the filaments in each direction is so thick that a motor 
can bind only to a filament of a given direction (or if the flexibility of the tail is limited). 
This coupled to inhomogeneities of the motor density is likely to lead to an x dependent 
N ± (see figure [1] where N ± is the integral of the motor density beneath the plus and 
minus filaments respectively). In reality it is probable that the situation is a mix with 
both a time and an x dependence with sensitivity to the details of both the structure 
of the bundles and the elasticity of the motor tail. In such evident below, the 
x dependence will dominate. 

When N ± (x) has an explicit x dependence the motion is more interesting. The 
bundle appears to be trapped for long times at certain locations (see figure [[2] for 
a sample trajectory obtained from numerics). When this happens the motion of the 
bundle appears oscillatory like even if there is no bistability (see figure [T2l green - C). 

Specifically N + (x) and N~(x) are now random variables with averages N + (x) and 
N~(x) respectively. The overline denotes an average over x locations (or equivalently 
an average over disorder realizations). In actual experiments it is more likely that a 
bundle composed of many randomly oriented filaments will not be symmetrical, namely 
N + (x) 7^ N~(x). Therefore N + (x) — N~(x), which signifies the overall polarity of the 
bundle, is expected from the central limit theorem to scale as y/£, where £ is overall 
length of the bundle. 

When N + ^ N~ , as discussed above, there is a non-zero mean velocity implying 
that iV + — N~ acts as an effective force. The random N + (x) — N~(x) thus induces a 
random effective force-field for the motion of the center of mass of the bundle. 

We expect for such problems, that on large-length scales and long-times, the systems 
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Figure 12: Representative, x(t), trajectories obtained from simulations with different 
values of fi. A is the distance between adjacent motors and a is a scale factor which 
allows us to show all three curves clearly on the same graph. In blue (A) a trajectory 
when there is a bistability, with /i ~ 0.4 and a = 1. In red (B) there is a bistability , 
with fi ~ 1 and a = 0.077. In green (C) there is no bistable region, with /i ~ 0.6 and 
a = 0.008. Inside boxes is a zoom on the motion. Parameters are the same as the ones 
in figures [TH and [T5j 



can be described by a particle diffusing in an effective random forcing field. This problem 
has been studied extensively in the past [28] and one finds four prominent types of 
behaviors: 

(1) fi > 2 ,(x(t)) ~ t and (x(t) 2 ) - (x(t)f ~ t. 

(2) 1 < fx < 2 - (x(t)} ~ t and (x(t) 2 ) - (x(t)f ~ t 2 ^. 

(3) < fi < 1 - (x(t)) ~ f and (x(t 2 )) - (x(t)} 2 ~ t 2 ^. 

(4) fi = - Sinai diffusion - F = 0, (x(t)) = and (a;(t) 2 ) ~ log^t/h). 

The exponent /x = 2F Q D/a is related to F , the average force acting on the particle, D, 
the diffusion constant in the absence of disorder and a is defined through the correlation 
F(x)F(x') = <j5(x — x'). The overline denotes an average over disorder realizations, 
angular brackets denote an average over histories of the system with a given realization 
of disorder and f(x) = F(x) + F is the force acting on the particle at x with the choice 
F — 0. Anomalous dynamics occur for fi < 2 and more prominently for fi < 1, and 
result from particles being trapped in rare deep wells, created by the random forcing 
field, for exponentially long times in the depth of the well [28] . 

To check if this straightforward analogy holds we account for disorder by defining a 
model, illustrated in figure [T3| where N ± (x) are random variables that depend on x, the 
location of the center of mass of the bundle. To generate a random motor landscape the 
motors are placed equally spaced, separated by the distance A, on a one dimensional 
lattice and assigned a random pulling direction that does not change with time. We 
consider a bundle composed of two groups of filaments with opposite polarities. The 
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Figure 13: The model we use to capture the effects of disorder. The bundle (yellow) 
is set on a surface of random motors with either a plus or minus pulling direction. In 
case (a) the pulling direction at each site changes randomly with time. In case (b) the 
pulling direction is chosen randomly and does not change with time. 



plus filament has a length of L + and the minus filament L_. Therefore the probability 
that a motor has a plus / minus pulling direction is p± = L±/£, where i = L + + L_. 
It is easy to verify that the average polarity is N + (x) — N~(x) oc L + — L~ . Note that 
since the motors are equally spaced the bundle is always subject to a constant number 
of motors N + (x) + N~(x) proportional to t. 

The scale of the effective random forcing energy landscape is set by the distance 
between adjacent motors A. The resulting effective potential is therefore expected to 
be proportional to A. In simulations /i can be adjusted by varying A while keeping all 
other parameters constant. This, however, may be difficult to achieve in experiments 
since A is controlled by the motor density. To keep all other parameters constant the 
size of the bundle and the viscosity must therefore also be scaled. Nonetheless the motor 
density is related to ji. 

We now turn to consider how a typical bundle will behave. A typical bundle 
composed of many filaments is expected to have an average polarity N + (x) — N~(x) oc 
VI. For this reason we set N + (x) — N (x) oc V N + + N~ . 

Using Monte-Carlo simulations and averaging over many realizations of disorder we 
verify that indeed the behavior is similar to the motion of a particle in a random forcing 
energy landscape and we are able to extract the exponent \i for different parameters 
from data over a few decades. We are interested in the dependence of fi on the bundle 
lengths £ in the bistable and non bistable regime. The results are given in figure UM and 
figure [15j We find two different behaviors for the bistable and non bistable regimes. For 
the bistable regime fi grows with iV + + N~ . Therefore in a motility assay experiment 
we expect that larger bundles will have larger [i (see figure (HJ). On long-time scales 
and large-length scales the motion will be a biased diffusion. On the other hand in the 
non bistable regime, as shown in figure [151 f 1 decreases as iV + + N~ increases. Here the 
x(t) grows slower for larger bundles and in the limit of an infinite bundle (x(t)) — > 
and (x 2 (t)) ~ log A (t/ti), with t\ a time scale. 

One characteristic of the motion in a random-forcing energy landscape with \i < 1 is 
the long period of time when the bundle is trapped in a "deep potential well" , localized 
around a certain point. These are then followed by a quick transition to another potential 
well (see figure fT2l) . On average the trapping time increases with the observation time. 
The farther the particle traverses the more likely it is to find a deeper well. Inside such 
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Figure 14: Disorder averaged trajectories, (x(t)), in the presence of a bistability. The 
black lines are the linear fits to the data which were used to obtain p. A = 0.25pm is 
the spacing between adjacent motors. In the inset p is shown as a function of N + + N~ . 
Here n/k a = 2, W = 10 (pm/sec) 2 , v = 2.5 pm/sec, G = 5 pN, r\ = 0.0125 pN sec/ pm 
and 7 = 1 pNsec/ pm. 
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Figure 15: Disorder averaged trajectories, (x(t)), when there is no bistability. The 
black lines are the linear fits to the data used to obtain p. A = 0.025/im is the spacing 
between adjacent motors. In the inset p is shown as a function of iV + + N~ and the 
line is a guide to the eye. n/k a = 2, W = 10 (pm/sec) 2 , v = 1 pm/sec, G = 5 pN, 
7] = 0.0125 pNsec/ pm and 7 = 1 pNsec/ pm. Parameters are the same as in figure [T41 
except for vq = 1 pm/sec. 
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a deep well the bundle is trapped and its motion is oscillatory like even if there is no 
bistability. Even for \x > 1 on short time scales the motion still includes long trapping 
times where the bundle is stuck at a certain point. 

We note that while in a tug-of-war where anomalous diffusion occurs naturally 
with no additional external force, it is possible to view the same effects by applying an 
external force to one set of motors. The force must then be tuned to reach the region 
of average zero velocity. 

6. Conclusions 

The paper focused on the dynamics of motors in a tug-of-war situation which are coupled 
by a rigid backbone. Using a steady-state mean-field solution of a model introduced 
in [20] we characterized different possible force-velocity relations which in many cases 
exhibit regions of bistability. The implications for motility assays were then discussed. 

It is interesting to compare our results to recent motility assays on myosin motors 
with a setup very similar to the one we consider [18]. In the experiments bundles of actin 
indeed exhibit trapping for long period of times in specific locations. When trapped they 
exhibit an oscillating like behavior between two velocities. This agrees well with our 
prediction for a random forcing energy landscape. 

Furthermore, in [18] they studied bundles trapped in such minima and considered 
the dependence of the average dwell time in each velocity as a function of the number 
of motors. The experiments show a weak dependence on the number of motors, which 
in [18] was explained by cooperative unbinding of the motors. As shown in figure [121 i n 
a minima of the random forcing energy landscape, as intuitively clear, the motion of the 
bundle exhibits an oscillating like behavior. Moreover, when the force- velocity relation 
exhibits no bistability the dwell time in each velocity has a very weak dependence 
on the number of motors (see figure [T6l) . This provides an alternative explanation of 
the experimental results. This conclusion could possibly be verified by single molecule 
experiments of the force-velocity curve of a single set of motors under the experimental 
settings of pS]. (We note that our analysis of figure [3] suggests a bistable behavior 
for myosin. However, as stated above, care has to be taken since a specific model was 
assumed in fitting the experimental data of [TT] and the exact setup in pj] and [18] 
is somewhat different.) For a clear analysis constant velocity experiments have to be 
carried out followed by the analysis we describe. In figure [TT] we show the average dwell 
time when a bistability is present in the force-velocity curve. There the behavior is 
distinct and the dwell time in each velocity has a strong exponential dependence on the 
number of motors. Note that the dwell time in figures [16] and [17| in principle depends 
on the realization of disorder and the local structure of the effective potential well. 

Finally, we note that using the techniques described above it is straight forward to 
derive the behavior of motors of different type each pulling in an opposite direction. 
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Figure 16: The dwell time, r, when there is no bistable region and N + — N~~ = 0. The 
dwell time is obtained from a single realization of the disorder using methods described 
in Appendix B As can be seen the dwell time has a weak dependence on N + + N~ . x(t) 
is sampled in steps of dt = 0.01/k. Here n/k a — 1, W — 10 (fim/sec) 2 , vq = 2 fim/sec, 
G = 5 pN, t] = 0.1 pNsec/ fim and 7 = 1.4 pN sec/ [im. 
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Figure 17: The dwell time, r, when there is a bistable region and N + — N~ = 0. The 
dwell time is obtained from a single realization of disorder using methods described 
in Appendix B As can be seen the dwell time, to a good approximation, grows 
exponentially with N + + N~ . The parameters are the same as in figure [16] except 
for 7 = 0.1 pNsec/ \im. 
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Appendix A. Simulations 

Monte Carlo simulations are employed to verify mean-field results and test our 
arguments. We use a standard Gillspie [29] algorithm. In each time step one transition 
occurs, namely one of N + + N~ motors will either detach or attach. Since transition 
times are Poissonian the probability for each transition is proportional to its respective 
rate. Next, the advancement time is chosen from an exponential distribution, with a 
time constant l/T = N+k d (v) + N~k d (-v) + (iV+ - N+)k a + (N~ - N~)k a . Once these 
two steps are completed the velocity and the rates are recalculated and the whole cycle 
is repeated. 

In the case of disorder the definition of the rates and the algorithm are mostly 
unchanged. The only difference is that one needs to account for the changes in N (x). 
This is done by modifying only the number of detached motors N^(x) as the bundle 
moves along the motor landscape. The results are not expected to change due to this 
approximation even under a more detailed model as long as most of the motors remain 
in the detached state. 

Throughout the paper x(t) is sampled at time steps of alt = 100/ac and the length 
of the simulations is 10 6 /k unless specified otherwise. 

Appendix B. Dwell times measurements 

In several parts of the paper the dwell times, r + and r_, are estimated in the bistable 
region from the numerical simulations. These are extracted from the autocorrelation 
function of N + — N~ which we observe to decay exponentially as, exp(— t/T) with T a 
time scale. In general r + ^ t_ and l/T — 1/t+ + 1/t_. To find r + and r_ independently 
we calculate numerically the ratio of the time spent in the positive and in the negative 
velocity which yields t + /t_. Then together with the expression for l/T r + and r_ are 
obtained. 

Note that when there is disorder the autocorrelation tail may not be purely 
exponential because of the local motor landscape. In this case the dwell time can 
be found directly from x(t). It is given by the average time it takes the velocity to 
switch signs. This method is useful when there is a single dwell time, namely f ext = 
and IF = IF. 
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